package column_generation;

import gurobi.GRB;
import gurobi.GRBEnv;
import gurobi.GRBException;
import gurobi.GRBLinExpr;
import gurobi.GRBModel;
import gurobi.GRBVar;

import Jama.Matrix;
import java.util.ArrayList;

public class RevisedSimplex {

    Matrix A;
    Matrix b;
    Matrix c;

    private final int b_length;
    private final int c_length;

    private int stockLength;
    private double[] specs;
    private boolean useGurobi = false;

    private Matrix B_inverse;
    private Matrix P;
    private Matrix E;
    private Matrix b_new;
    private ArrayList<Integer> bv;
    private boolean isOptimal = false;

    RevisedSimplex(Matrix A, Matrix b, Matrix c) {
        this.A = A;
        this.b = b;
        this.c = c;
        this.P = A.copy();
        this.b_length = b.getRowDimension();
        this.c_length = c.getColumnDimension();
    }

    RevisedSimplex(Matrix A, Matrix b, Matrix c, int stockLength, double[] specs) {
        this(A, b, c);
        this.stockLength = stockLength;
        this.specs = specs;
        this.useGurobi = true;
    }

    public void showModel() {
        // 输出目标函数
        String obj = "Min z = ";
        for (int i = 0; i < c_length; i++) {
            obj += c.get(0, i) + "x" + Integer.toString(i);
            if (i != c_length - 1) {
                obj += " + ";
            }
        }
        System.out.println(obj);
        System.out.println("Subject to:");

        // 输出约束
        for (int i = 0; i < b_length; i++) {
            String con = "";
            for (int j = 0; j < c_length; j++) {
                con += Double.toString(A.get(i, j)) + "x" + Integer.toString(j);
                if (j != c_length - 1) {
                    con += " + ";
                }
            }
            con += " = " + Double.toString(b.get(i, 0));
            System.out.println(con);
        }
    }

    private Matrix getCbv() {
        Matrix temp = new Matrix(1, b_length);
        for (int i = 0; i < b_length; i++) {
            temp.set(0, i, c.get(0, i));
        }
        return temp;
    }

    private int getXin() {
        int xin = 0;

        // 计算检验数
        double cmin = Double.MAX_VALUE;
        for (int i = 0; i < c_length; i++) {
            // 计算检验数

            double temp = c.get(0, i) - getCbv().times(P.getMatrix(0, b_length - 1, i, i)).det();
            if (cmin > temp) {
                cmin = temp;
                xin = i;
            }
        }
        if (cmin >= 0) {
            isOptimal = true;
        }
        return xin;
    }

    private int getXinByGurobi() {
        int xin = 0;
        Matrix p = new Matrix(b_length, 1);

        try {

            // Create empty environment, set options, and start
            GRBEnv env = new GRBEnv(true);
            env.start();

            // Create empty model
            GRBModel model = new GRBModel(env);

            // Create variables
            GRBVar[] x = model.addVars(b_length, GRB.INTEGER);

            // Set objective: maximize x + y + 2 z
            GRBLinExpr expr = new GRBLinExpr();
            expr.addTerms(getCbv().times(B_inverse).getArray()[0], x);
            expr.addConstant(-1);
            model.setObjective(expr, GRB.MAXIMIZE);

            // Add constraint
            expr = new GRBLinExpr();
            expr.addTerms(specs, x);
            model.addConstr(expr, GRB.LESS_EQUAL, stockLength, "c0");

            // Optimize model
            model.optimize();

            for (int i = 0; i < b_length; i++) {
                p.set(i, 0, x[i].get(GRB.DoubleAttr.X));
            }
            if (model.get(GRB.DoubleAttr.ObjVal) <= 0) {
                isOptimal = true;
            }

            // Dispose of model and environment
            model.dispose();
            env.dispose();

        } catch (GRBException e) {
            System.out.println("Error code: " + e.getErrorCode() + ". " + e.getMessage());
        }
        for (int i = 0; i < c_length; i++) {
            if (A.getMatrix(0, b_length - 1, i, i).minus(p).norm1() == 0) {
                xin = i;
            }
            // boolean flag = false;
            // for (int j = 0; j < b_length; j++) {
            // if (A.get(j, i) != p.get(j, 0)) {
            // flag = true;
            // break;
            // }
            // }
            // if (!flag) {
            // xin = i;
            // }
        }
        return xin;
    }

    private int getXout(int xin) {
        int xout = 0;
        double theta = Double.MAX_VALUE;
        // 计算theta
        for (int i = 0; i < b_length; i++) {
            double temp = b_new.get(i, 0) / P.getMatrix(0, b_length - 1, xin, xin).get(i, 0);
            if (temp < 0) {
                continue;
            }
            if (temp < theta) {
                theta = temp;
                xout = bv.get(i);
            }
        }
        return xout;
    }

    private void updateB_inverse(int xin, int xout) {
        Matrix p = P.getMatrix(0, b_length - 1, xin, xin);
        // 构造E
        int index = bv.indexOf(xout);
        for (int i = 0; i < b_length; i++) {
            Matrix temp = new Matrix(b_length, 1);
            if (i == index) {
                for (int j = 0; j < b_length; j++) {
                    if (j == index) {
                        temp.set(j, 0, 1 / p.get(j, 0));
                    } else {
                        temp.set(j, 0, -p.get(j, 0) / p.get(index, 0));
                    }
                }
            } else {
                for (int j = 0; j < b_length; j++) {
                    if (j == i) {
                        temp.set(j, 0, 1);
                    } else {
                        temp.set(j, 0, 0);
                    }
                }
            }
            E.setMatrix(0, b_length - 1, i, i, temp);
        }
        // 更新B_inverse
        B_inverse = E.times(B_inverse);
        b_new = B_inverse.times(b);
    }

    void solve() {
        // 确定初始基变量
        bv = new ArrayList<Integer>(b_length);
        for (int i = 0; i < b_length; i++) {
            bv.add(-1);
        }
        for (int i = 0; i < c_length; i++) {
            Matrix pTemp = A.getMatrix(0, b_length - 1, i, i);
            int[] temp = new int[b_length];
            int p = -1; // 指向某一列最后一个不为0的元素
            for (int j = 0; j < b_length; j++) {
                if (pTemp.get(j, 0) != 0) {
                    temp[j] = 1;
                    p = j;
                }
            }
            int sum = 0;
            for (int j = 0; j < b_length; j++) {
                sum += temp[j];
            }
            if (sum == 1) {
                bv.set(p, i);
            }
        }

        // （假装）计算初始B逆
        int[] bvList = new int[b_length];
        for (int i = 0; i < b_length; i++) {
            bvList[i] = bv.get(i);
        }
        B_inverse = A.getMatrix(0, b_length - 1, bvList).inverse();
        b_new = B_inverse.times(b);

        // 生成初始E（单位矩阵）
        E = new Matrix(b_length, b_length);
        for (int i = 0; i < b_length; i++) {
            for (int j = 0; j < b_length; j++) {
                if (i == j) {
                    E.set(i, j, 1);
                } else {
                    E.set(i, j, 0);
                }
            }
        }

        // 开始计算
        long startTime = System.currentTimeMillis();
        while (true) {
            B_inverse.print(3, 3);
            System.out.println("基变量为：" + bv.toString());
            // 更新变量技术参数向量
            for (int i = 0; i < c_length; i++) {
                Matrix p = new Matrix(b_length, 1);
                for (int j = 0; j < b_length; j++) {
                    p.set(j, 0, A.get(j, i));
                }
                P.setMatrix(0, b_length - 1, i, i, B_inverse.times(p));
            }
            // 开始计算
            int xin = 0;
            if (useGurobi) {
                xin = getXinByGurobi();
                System.out.println(xin);
            } else {
                xin = getXin();
            }
            // 输出当前解
            System.out.println("Z = " + Double.toString(getCbv().times(b_new).det()));
            // 判断是否得到最优解
            if (isOptimal) {
                System.out.println("最优解为：");
                for (int i = 0; i < c_length; i++) {
                    System.out.println("x" + Integer.toString(i) + " = "
                            + (bv.indexOf(i) != -1 ? Double.toString(b_new.get(bv.indexOf(i), 0)) : "0"));
                }
                break;
            }
            int xout = getXout(xin);
            updateB_inverse(xin, xout);
            bv.set(bv.indexOf(xout), xin);
            System.out.println("==========================");
        }
        long endTime = System.currentTimeMillis();
        System.out.println("程序运行时间：" + (endTime - startTime) + "ms");
    }

}
